ptm = proc.time()

1 Introdução

Frequentemente utilizamos serviços de recomendação seja por plataformas como Netflix, Mercado Livre, OLX dentre outros. O conjunto de dados explanados a seguir foram extraídos da plataforma Kaggle e são referentes ao catálogo de filmes da plataforma IMDB, o intuito aqui é criar um sistema de classificação de estrelas orientado a dados do catálogo e avaliação de filmes do IMDB.

Esse projeto irá explorar o conjunto de dados e gerar insights, e aplicar o algoritmo de Random Forest (baseado na clusterização de Breiman).

Esse trabalho irá levantar quais os fatores mais importantes para que um filme tenha uma alta classificação (categoria de faixas de estrelas), utilizando um modelo de Random Forest para classificação e resultados escritos em R.

Pacotes e Leitura

Alguns pacotes R utilizados nessa rotina

library(tidyverse)          # Manipulação de banco de dados e análise exploratória
library(highcharter)        # Highchart gráficos usando htmlwidgets (renderiza HTML via código R Markdown)
library(DT)                 # Tabelas: datatable()
library(knitr)          
library(kableExtra)
library(rpart.plot)
library(randomForest)
library(corrplot)
library(wordcloud2)
RECOMEND.METADATA = readxl::read_xlsx(path = ".../DB/IMDB.xlsx", sheet = 1)
db = RECOMEND.METADATA

2 Banco de Dados

O conjunto de dados utilizado foi o de título IMDB 5000 extraído da plataforma kaggle. As informações contidas no banco foram catalogadas de filmes publicados ao longo de 100 anos em 66 países (entre 1916 e 2016) da plataforma IMDB - Internet Movie DataBases, o Arquivo original contém 5044 filmes (observações) e 28 variáveis descritas a seguir.

Os dados em questão são públicos e disponíveis para download clicando AQUI ou AQUI.

3 Análise Exploratória dos dados

## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

3.1 Estrutura dos dados

## O banco de dados possui 5043 observações e 28 variáveis

Excluindo algumas variáveis que não serão utilizadas.

    1. As variáveis color, director_name, language, country, content_rating, em especial, foram removidas uma vez que estas não se mostraram importantes no algoritmo de importância do Random Forest de Breiman.
    1. Já as variáveis director_facebook_likes, actor_3_facebook_likes, actor_2_name, actor_1_facebook_likes, actor_1_name, actor_3_name,-facenumber_in_poster, actor_2_facebook_likes, movie_imdb_link foram extraídas de forma determinística por escolha dos autores deste.

O motivo de remover as variáveis em (1) é devido à tentativa de reduzir o número de missing na base de dados. no passo 3.1 de Tratamento de NA’s.

#db=RECOMMEND.METADATA
db = db %>% select(- director_facebook_likes, -actor_3_facebook_likes, -actor_2_name, -actor_1_facebook_likes, -actor_1_name, -actor_3_name,-facenumber_in_poster, -actor_2_facebook_likes, -color, -director_name, -language, -country, -content_rating, -movie_imdb_link)
## O banco de dados agora possui 14 variáveis. Ou seja, foram removidas  14 variáveis no passo anterior.

Estrutura dos dados

glimpse(db) #str(db)
## Observations: 5,043
## Variables: 14
## $ num_critic_for_reviews    <int> 723, 302, 602, 813, NA, 462, 392, 32...
## $ duration                  <int> 178, 169, 148, 164, NA, 132, 156, 10...
## $ gross                     <int> 760505847, 309404152, 200074175, 448...
## $ genres                    <fct> Action|Adventure|Fantasy|Sci-Fi, Act...
## $ movie_title               <fct> Avatar , Pirates of the Caribbean: A...
## $ num_voted_users           <int> 886204, 471220, 275868, 1144337, 8, ...
## $ cast_total_facebook_likes <int> 4834, 48350, 11700, 106759, 143, 187...
## $ plot_keywords             <fct> avatar|future|marine|native|parapleg...
## $ num_user_for_reviews      <int> 3054, 1238, 994, 2701, NA, 738, 1902...
## $ budget                    <dbl> 237000000, 300000000, 245000000, 250...
## $ title_year                <int> 2009, 2007, 2015, 2012, NA, 2012, 20...
## $ imdb_score                <fct> 7.9, 7.1, 6.8, 8.5, 7.1, 6.6, 6.2, 7...
## $ aspect_ratio              <fct> 1.78, 2.35, 2.35, 2.35, , 2.35, 2.35...
## $ movie_facebook_likes      <int> 33000, 0, 85000, 164000, 0, 24000, 0...

3.2 Transformação de variáveis

Transformando as variáveis imdb_score e aspect_ratio em numéricas

db$imdb_score = as.numeric(as.character(db$imdb_score))
db$aspect_ratio = as.numeric(as.character(db$aspect_ratio))
hc1 = hchart(db$imdb_score, color = "#e8bb0b", name = "imdb_score") %>% 
        hc_title(text = "Histograma dos votos na plataforma IMDB") %>%
        hc_exporting(enabled = TRUE, filename = "Fig1-Pimenta"); hc1
hc2 = hchart(db$num_voted_users, color = "#786eea", name = "imdb_num_votos") %>% 
        hc_title(text = "Histograma dos número de votos na plataforma IMDB") %>%
        hc_exporting(enabled = TRUE, filename = "Fig2-Pimenta"); hc2

Criando uma função para cálculo da moda

Moda <- function(x) {
     ux <- unique(x)
     ux[which.max(tabulate(match(x, ux)))]
}
## Média, moda e mediana da variável imdb_score são, respectivamente, 6.44 6.7 6.6

Extraindo valores da variável gênero e transformando em dummies

gg <- as.character(db$genres); gg <- gsub("-", "_", as.character(gg))

#t <- unlist(strsplit(gg[1],split = "\\|"))
tem1 <- data.frame() 
for(i in 1:length(gg)){
        tem <- tem1
        t <- unlist(strsplit(gg[i],split = "\\|"))
        temp <- data.frame(t)
        tem1 <- rbind(tem,temp)
}

Gen <- unique(tem1); Gen <- gsub("-", "_", as.character(Gen$t))
cat("Existem", length(Gen), "valores de gêneros únicos de filme no banco de dados.")
## Existem 26 valores de gêneros únicos de filme no banco de dados.
Genname <- Gen

fe <- matrix(data = 0, nrow = length(gg), ncol = length(Genname))
fe <- data.frame(fe); colnames(fe) <- Genname

for(i in 1:length(gg)){
        for(j in 1:length(Genname)){
                g <- grepl(Genname[j], gg[i])
                if(g == TRUE){
                        fe[i, j] <- 1        
                }
        }
}

NumGen = as_tibble(rbind(apply(fe,2,sum)))
NumGen = gather(NumGen, key = "variables", value = "num_gender")
NumGen = NumGen[order(NumGen$num_gender, decreasing = TRUE), ]


hc3 <- highchart() %>%
  hc_add_series(data = NumGen$num_gender, 
                type = "bar",
                name = "# de filmes",
                showInLegend = FALSE,
                tooltip = list(valueDecimals = 0, valuePrefix = "", valueSuffix = ""), color="blue") %>%
  hc_yAxis(title = list(text = "Quantitativo de filmes"), 
           allowDecimals = TRUE, max = (max(NumGen$num_gender)+103),
           labels = list(format = "{value}")) %>%
  hc_xAxis(title = list(text = "Gênero de filme"),
           categories = NumGen$variables,
           tickmarkPlacement = "on",
           opposite = FALSE) %>%
  hc_title(text = "Quantitativo de filmes por gênero",
           style = list(fontWeight = "bold")) %>% 
  hc_subtitle(text = paste("")) %>%
      hc_tooltip(valueDecimals = 2,
                 pointFormat = "{point.y} filmes")%>%
                 #pointFormat = "Variável: {point.x} <br> Missing: {point.y}") 
      hc_credits(enabled = TRUE, 
                 text = "Fonte: IMDB/KAGGLE. Elaboração: Ewerson Pimenta.",
                 style = list(fontSize = "10px")) %>%
  hc_exporting(enabled = TRUE, filename = "F3-filmes-genero-Pimenta")
#hc <- hc %>% 
#  hc_add_theme(hc_theme_darkunica())
hc3; remove(gg,t,tem, temp, tem1,Gen, g, i, j)

Unificando a base de dados com as novas variáveis de gêneros únicos descobertos.

db1 = as_tibble(cbind(as.data.frame(db),fe))
cat("Foram inseridas todas as", dim(fe)[2],"novas variáveis provenientes dos gêneros únicos descobertos no passo anterior.")
## Foram inseridas todas as 26 novas variáveis provenientes dos gêneros únicos descobertos no passo anterior.

Algumas variáveis de gêneros possuem pouquíssimos filmes e também serão removidas do banco. São elas Film_Noir, Short, News, Reality_TV e Game_Show.

db1$Film_Noir <- NULL
db1$Short <- NULL
db1$News <- NULL
db1$Reality_TV <- NULL
db1$Game_Show <- NULL

Semelhante à variável de gênero, foi feito o split por palavra chave da variável keyword, esse, por sua vez, continha 6780 palavras únicas. E então, geramos uma nuvem de palavras chave, com o auxílio do pacote wordcloud2 com o seguinte comando, após obter o data frame do nome das palavras chaves e frequência.

library(wordcloud2)
wordcloud2(NumKeyWord)

3.3 Tratamento de NA e duplicatas

Porcentagem de NA por variável

db_miss <- db %>% summarise_all(funs(sum(is.na(.))/n()))
db_miss <- gather(db_miss, key = "variables", value = "percent_missing")
db_miss$percent_missing = 100*db_miss$percent_missing
db_miss = db_miss[order(db_miss$percent_missing, decreasing = TRUE), ]
#db_miss

hc4 <- highchart() %>%
  hc_add_series(data = db_miss$percent_missing, 
                type = "bar",
                name = "Porcentagem de missing",
                showInLegend = FALSE,
                tooltip = list(valueDecimals = 2, valuePrefix = "", valueSuffix = " %")) %>%
  hc_yAxis(title = list(text = "Porcentagem de missing"), 
           allowDecimals = TRUE, max = 20,
           labels = list(format = "{value}%")) %>%
  hc_xAxis(categories = db_miss$variables,
           tickmarkPlacement = "on",
           opposite = FALSE) %>%
  hc_title(text = "Porcentagem de missinig por variável",
           style = list(fontWeight = "bold")) %>% 
  hc_subtitle(text = paste("")) %>%
      hc_tooltip(valueDecimals = 2,
                 pointFormat = "Missing: {point.y}")%>%
                 #pointFormat = "Variável: {point.x} <br> Missing: {point.y}") 
      hc_credits(enabled = TRUE, 
                 text = "Fonte: IMDB/KAGGLE. Elaboração: Ewerson Pimenta.",
                 style = list(fontSize = "10px")) %>%
  hc_exporting(enabled = TRUE, filename = "Fig0-Pimenta")
#hc <- hc %>% 
#  hc_add_theme(hc_theme_darkunica())
hc4

Removendo todas as observações que contêm NA.

db1 <- db1 %>% drop_na(gross, budget, aspect_ratio 
,title_year 
,num_critic_for_reviews, num_user_for_reviews)
#glimpse(db1)

db1 <- db1 %>% na.omit()
db1 <- db1 %>% drop_na()

Removendo dados duplicados

cat("Existem", sum(duplicated(db1)), "filmes duplicados na base de dados.")
## Existem 33 filmes duplicados na base de dados.
db1 <- db1[!duplicated(db1), ]

Excluindo filmes lançados antes de 1980

hc00 = hchart(db1$title_year, color = "#53a074", name = "imdb1_score") %>% 
        hc_title(text = "Histograma do ano de publicação dos filmes") %>%
        hc_exporting(enabled = TRUE, filename = "Fig1-Pimenta"); hc00

Percebemos pelo gráfico anterior que há muitos poucos filmes publicados antes de 1980 e estes podem não ser representativos. Decidimos, então, por trabalhar apenas com os filmes que foram publicados a aprtir de 1980

db1 <- db1[db1$title_year >= 1980,]
hc01 = hchart(db1$title_year, color = "#79d8a2", name = "imdb1_score") %>% 
        hc_title(text = "Histograma do ano de publicação dos filmes") %>%
        hc_exporting(enabled = TRUE, filename = "Fig1-Pimenta"); hc01
## O novo banco de dados, sem observações faltantes, possui 3687 observações. Ou seja, no processo de remoção de valores faltantes foram perdidas 1356 observações. Finalmente, removemos todas as observações duplicadas e faltantes do banco.

3.4 Weighted Rating (WR) - IMDB_score

Um passo importante é penalizar a variável de escore imdb_score pelo número de votos recebidos. Ver mais no estimador de Shrinkage

\(WR = \frac{v}{v+m} \times R + \frac{m}{v+m} \times C\)

Onde,

R = as.numeric(db1$imdb_score)
v = as.numeric(db1$num_voted_users)
m = 7000
C = mean(db1$imdb_score)
db1$WR = (v/(v+m))*R + (m/(v+m))*C
  • R = Escore médio dos votos para o título do filme dado pelos usuários do IMDB = (imdb_score)
  • v = Número de usuários que votaram = (num_voted_users)
  • m = Mínimo de votos requerido (atualmente 7.000)
  • C = O escore médio de todos os 3766 filmes (atualmente 6,5)

Medidas pontuais e de dispersão de imdb_score e WR

summary(db1$imdb_score)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.600   5.900   6.500   6.442   7.200   9.300
summary(db1$WR)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.017   6.031   6.496   6.514   7.053   9.288
cat("Além disso, conseguimos reduzir o desvio padrão que orignalmente era de ",round(sd(db1$imdb_score),2),",para a variável imdb_score, e agora passou a ser ",round(sd(db1$WR),2), "para a variável WR.")
## Além disso, conseguimos reduzir o desvio padrão que orignalmente era de  1.05 ,para a variável imdb_score, e agora passou a ser  0.84 para a variável WR.

Selecionando uma amostra aleatória (a.a.) de tamanho \(n=800\) e representando em um gráfico de dispersão com valores reais vs ajustados.

set.seed(123654)
amostra0 = sample(x = 1:dim(db1)[1], size = 800, replace = FALSE)
dbX = db1[amostra0,] %>% 
  select(movie_title,num_voted_users,imdb_score, WR)
  
dss <- map(c("cross"), function(s){
  
  x <- as.numeric(dbX$imdb_score)
  y <- as.numeric(dbX$WR)
  
  list(name = s,
       data = list_parse(data_frame(x, y)),
       marker = list(symbol = s, enabled = TRUE), lineColor = "#56667a")
  
})
#dss[[1]]$data[amostra1]

hc5 = highchart() %>% 
  hc_chart(type = "scatter", color = "#56667a") %>% 
  hc_title(text = "Score IMDB vs WR (calculado pelo estimador de Shrinkage)") %>%
  hc_subtitle(text = "800 filmes selecionados via amostra aleatória simples") %>%
  hc_xAxis(title = list(text = "x: imdb_score"), 
           allowDecimals = TRUE, labels = list(format = "{value}★")) %>%
  hc_yAxis(title = list(text = "y: WR (calibrado)"),
           allowDecimals = TRUE, labels = list(format = "{value}★")) %>%
  hc_exporting(enabled = TRUE, filename = "F3-Pimenta") %>%
  hc_add_series_list(dss); hc5; remove(dbX, dss)

Análise do cálculo de IMDB

  • Para filmes com # de votos recebidos MENOR que 7mil (m: votos mínimos requeridos):
    • imdb_score > 6,5: DECRESCIMENTO
    • imdb_score < 6,5: CRESCIMENTO
  • Para filmes com # de votos recebidos MAIOR que 7mil (m: votos mínimos requeridos):
    • imdb_score > 6,5: DECRESCIMENTO
    • imdb_score < 6,5: CRESCIMENTO

Ou seja, para os filmes catalogados com m muito inferior a 7mil o novo escore calibrado teve maior diferentça que aqueles superior a 7 mil. Além disso, quando scores são maiores que 6,5 o WR tende a cair, caso contrário o valor pode aumentar. Quanto maior o número de votos recebidos, menor a diferença do valor de IMDB_score e WR.

set.seed(123654)
amostra2 = sample(x = 1:dim(db1)[1], size = 35, replace = FALSE)
db1[amostra2,] %>% 
  select(movie_title,num_voted_users,imdb_score, WR, movie_facebook_likes, budget) %>% 
  datatable(rownames = amostra2,options = list(searchin = TRUE, scrollX = TRUE, pageLength = 5))

3.5 Visualizações Finais

Finalmente, segue uma a.a.da base de dados após limpeza e análise exploratória.

set.seed(123851)
amostra3 = sample(x = 1:dim(db1)[1], size = 35, replace = FALSE)
db1[amostra3,] %>%
datatable(rownames = amostra3,options = list(searchin = TRUE, scrollX = TRUE, pageLength = 5))

Para os dados quantitativos a seguinte Matriz de correlação. Foi gerada utilizando o pacote corrplot

#install.packages("corrplot")
library(corrplot)  
res1 <- cor.mtest(db_cor, conf.level = .95)
res2 <- cor.mtest(db_cor, conf.level = .99)

corrplot::corrplot(cor(db_cor), method = "color", col = col(200),
         type = "upper", order = "hclust", number.cex = .7,
         addCoef.col = "black", # Add coefficient of correlation
         tl.col = "black", tl.srt = 90, # Text label color and rotation
         # Combine with significance
         p.mat = res1$p, sig.level = 0.01, insig = "blank", 
         # hide correlation coefficient on the principal diagonal
         diag = FALSE)

4 Modelagem

4.1 Preparação de dados

Transformação da variável contínua WRem categórica WR_Grp

Transformando a variável WR em fator, com as seguintes categorias:

  • 0 a 4 ★
  • de 4 a 6 ★
  • de 6 a 8 ★
  • de 8 a 10 ★
movie = db1
Grp <- function(tn){
  tn = abs(tn)
    if (tn >= 0 & tn <= 4){
        return('de 0 a 4')
    }else if(tn > 4 & tn <= 6){
        return('de 4 a 6')
    }else if(tn > 6 & tn <= 8){
        return('de 6 a 8')
    }else if (tn > 8){
        return('de 8 a 10')
    }
}
# apply the Group function to the WR column
movie$WR_Grp <- sapply(movie$WR,Grp)
# set as factor the new column
movie$WR_Grp <- as.factor(movie$WR_Grp)
#View(head(movie))
table(movie$WR_Grp)
## 
##  de 0 a 4  de 4 a 6  de 6 a 8 de 8 a 10 
##        29       843      2691       124
# apply the Group function to the WR column
imdb_score_Grp <- sapply(movie$imdb_score,Grp)
# set as factor the new column
imdb_score_Grp <- as.factor(imdb_score_Grp)
#View(head(movie))
table(imdb_score_Grp)
## imdb_score_Grp
##  de 0 a 4  de 4 a 6  de 6 a 8 de 8 a 10 
##        92      1049      2413       133

Remoção de variáveis Remove as variáveis imdb_score, genres, plot_keywords, movie_imdb_link e WR.

movie$imdb_score <- NULL
movie$genres <- NULL
movie$plot_keywords <- NULL
movie$movie_imdb_link <- NULL
movie$WR <- NULL
glimpse(movie)
## Observations: 3,687
## Variables: 33
## $ num_critic_for_reviews    <int> 723, 302, 602, 813, 462, 392, 324, 6...
## $ duration                  <int> 178, 169, 148, 164, 132, 156, 100, 1...
## $ gross                     <int> 760505847, 309404152, 200074175, 448...
## $ movie_title               <fct> Avatar , Pirates of the Caribbean: A...
## $ num_voted_users           <int> 886204, 471220, 275868, 1144337, 212...
## $ cast_total_facebook_likes <int> 4834, 48350, 11700, 106759, 1873, 46...
## $ num_user_for_reviews      <int> 3054, 1238, 994, 2701, 738, 1902, 38...
## $ budget                    <dbl> 237000000, 300000000, 245000000, 250...
## $ title_year                <int> 2009, 2007, 2015, 2012, 2012, 2007, ...
## $ aspect_ratio              <dbl> 1.78, 2.35, 2.35, 2.35, 2.35, 2.35, ...
## $ movie_facebook_likes      <int> 33000, 0, 85000, 164000, 24000, 0, 2...
## $ Action                    <dbl> 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, ...
## $ Adventure                 <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ Fantasy                   <dbl> 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, ...
## $ Sci_Fi                    <dbl> 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, ...
## $ Thriller                  <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Documentary               <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Romance                   <dbl> 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, ...
## $ Animation                 <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, ...
## $ Comedy                    <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, ...
## $ Family                    <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, ...
## $ Musical                   <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, ...
## $ Mystery                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ...
## $ Western                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Drama                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ History                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Sport                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Crime                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Horror                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ War                       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Biography                 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Music                     <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, ...
## $ WR_Grp                    <fct> de 6 a 8, de 6 a 8, de 6 a 8, de 8 a...
cat("A base de dados que iremos trabalhar no modelo tem", dim(movie)[1] , "observações e ", dim(movie)[2]-1, "variáveis, sem contar a variável que contém os títulos dos filmes.")
## A base de dados que iremos trabalhar no modelo tem 3687 observações e  32 variáveis, sem contar a variável que contém os títulos dos filmes.

Partição dos dados

Particionando a base de dados em Treino e Teste, esses dois (Treino e Teste) também terão armazenos os nomes dos filmes selecionados via amostragem probabilística dos dados originais separadamente das bases de Treino e Teste.

Posteriormente, removemos os rótulos dos filmes nas bases Treino e Teste

set.seed(9182345)
ind <- sample(2, nrow(movie), replace = T, prob = c(0.7, 0.3))
train <- movie[ind==1, -4]
test <- movie[ind==2, -4]
trainMovie <- movie[ind==1, 4]
testMovie <- movie[ind==2, 4]

Distribuição dos escores nas bases de treino e teste

round(table(train$WR_Grp)/sum(table(train$WR_Grp))*100,2)
## 
##  de 0 a 4  de 4 a 6  de 6 a 8 de 8 a 10 
##      0.88     22.70     72.69      3.73
round(table(test$WR_Grp)/sum(table(test$WR_Grp))*100,2)
## 
##  de 0 a 4  de 4 a 6  de 6 a 8 de 8 a 10 
##      0.55     23.25     73.71      2.49

4.2 Random Forest - Metodologia

Descrição 1. Random Forest foi desenvolvido para agregar árvores de decisão (modelo de classificação);
2. Pode ser usado para modelo de classificação (p/ var. resposta categórica) ou regressão (no caso de haver variável resposta contínua);
3. Evita overfitting;
4. Permite trabalhar com um largo número de características de um conjunto de dados;
5. Auxilia na seleção de variáveis baseada em um algoritmo que calcula a importância por variável (assim, tendo conhecimento de quais variáveis são mais importantes, podemos usar essa informação para outros modelos de classificação);
6. User-friendly: apenas 2 parâmetros livres:

  • Trees - ntrees, default 500 (Nº de árvores);
  • Variáveis selecionadas via amostragem aleatória candidatas à cada “split” (quebra da árvore) - mtry, default \(\sqrt{p}\) p/ classificação e \(\frac{p}{3}\) p/ regressão (p: nº de features/variáveis);

Passo-a-Passo

É realizado em 3 passos:

  1. Desenha as amostras via bootstrap do número de árvores ntrees;
  2. Para cada amostra via bootstrap, cresce o número de árvores “un-puned” para a escolha da melhor quebra da árvore baseado na amostra aleatória do valor predito de mtry a cada nó da árvore;
    1. Faz classificação de novos valores usando a maioria de votos p/ classificação e usa a média p/ regressão baseada nas amostras de ntrees.

Exemplo

4.3 Random Forest - Aplicação e Resultados

Inicialmente utilizaremos o pacote randomForest que implmenta o algoritmo de Random Forest de Breiman (baseado na clusterização de Breiman, originalmente codificada em Fortran) que tem por finalidade classificar e/ou criar regressão. Além disso, pode ser usado em um modelo não supervisionado para avaliar proximidades entre pontos.

Estamos usando, a partir daqui, a base de treino.

#library(randomForest)
#library(rpart)
#library(rpart.plot)
#rf <- randomForest(proximity = T,ntree = 38,do.trace = T,WR~.,data=training)
set.seed(9984512)
rf <- randomForest(WR_Grp~.,data=train)
rf
## 
## Call:
##  randomForest(formula = WR_Grp ~ ., data = train) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 5
## 
##         OOB estimate of  error rate: 18.02%
## Confusion matrix:
##           de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10 class.error
## de 0 a 4         0       11       12         0  1.00000000
## de 4 a 6         0      291      300         0  0.50761421
## de 6 a 8         0       95     1790         7  0.05391121
## de 8 a 10        0        0       44        53  0.45360825
attributes(rf)
## $names
##  [1] "call"            "type"            "predicted"      
##  [4] "err.rate"        "confusion"       "votes"          
##  [7] "oob.times"       "classes"         "importance"     
## [10] "importanceSD"    "localImportance" "proximity"      
## [13] "ntree"           "mtry"            "forest"         
## [16] "y"               "test"            "inbag"          
## [19] "terms"          
## 
## $class
## [1] "randomForest.formula" "randomForest"

Olhando as 6 primeiras observações real X predito

p1 <- predict(rf,train)
head(p1)
##         1         2         3         4         5         6 
##  de 6 a 8  de 6 a 8  de 6 a 8 de 8 a 10  de 6 a 8  de 6 a 8 
## Levels: de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10
head(train$WR_Grp)
## [1] de 6 a 8  de 6 a 8  de 6 a 8  de 8 a 10 de 6 a 8  de 6 a 8 
## Levels: de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10

Matriz de confusão

library(caret)
library(e1071)
confusionMatrix(p1, train$WR_Grp)
## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10
##   de 0 a 4        23        0        0         0
##   de 4 a 6         0      591        0         0
##   de 6 a 8         0        0     1892         0
##   de 8 a 10        0        0        0        97
## 
## Overall Statistics
##                                      
##                Accuracy : 1          
##                  95% CI : (0.9986, 1)
##     No Information Rate : 0.7269     
##     P-Value [Acc > NIR] : < 2.2e-16  
##                                      
##                   Kappa : 1          
##  Mcnemar's Test P-Value : NA         
## 
## Statistics by Class:
## 
##                      Class: de 0 a 4 Class: de 4 a 6 Class: de 6 a 8
## Sensitivity                 1.000000           1.000          1.0000
## Specificity                 1.000000           1.000          1.0000
## Pos Pred Value              1.000000           1.000          1.0000
## Neg Pred Value              1.000000           1.000          1.0000
## Prevalence                  0.008836           0.227          0.7269
## Detection Rate              0.008836           0.227          0.7269
## Detection Prevalence        0.008836           0.227          0.7269
## Balanced Accuracy           1.000000           1.000          1.0000
##                      Class: de 8 a 10
## Sensitivity                   1.00000
## Specificity                   1.00000
## Pos Pred Value                1.00000
## Neg Pred Value                1.00000
## Prevalence                    0.03726
## Detection Rate                0.03726
## Detection Prevalence          0.03726
## Balanced Accuracy             1.00000

RF_importance = randomForest::importance(rf)[order(randomForest::importance(rf)[,1], decreasing = TRUE), ]
knitr::kable(RF_importance)
x
num_voted_users 165.997906
duration 101.645711
gross 92.466271
num_user_for_reviews 89.211429
num_critic_for_reviews 85.536687
budget 85.152515
cast_total_facebook_likes 77.366530
title_year 76.802149
movie_facebook_likes 63.545716
Drama 46.647851
Horror 24.058601
aspect_ratio 16.664637
Action 15.235928
Comedy 15.109218
Thriller 11.041651
Sci_Fi 10.834496
Fantasy 10.562471
Romance 9.785667
Adventure 9.376240
Crime 9.300023
Family 9.230662
Animation 8.317384
Mystery 6.955240
Music 5.951651
Biography 3.720929
Sport 3.365780
Musical 2.598940
War 2.515458
Western 1.767951
Documentary 1.466404
History 1.426391

Verificamos que as variávei Game_Show, Sci_Fi, Reality_TV, News e Film_Noir não foram relevantes para o algoritmo do random forest.

randomForest::varImpPlot(rf)

Taxa de Erro - Random Forest

plot(rf)
legend('topright', colnames(rf$err.rate), col=1:5, fill=1:5)

Observamos que, a partir do número de árvores geradas \(ntrees > 100\) o erro OOB (Out of Bag) não pode ser melhorado.

Ajustando e melhorando estimativas

Além disso, iremos alterar alguns parâmetros da função randomForest como o número de ntrees e mtry. Assim, repetimos o algoritmo do Random Forest, ainda usando a base treino.

# Tune mtry
x = as.data.frame(train1[,-31])
y = (as.factor(train1$WR_Grp))
t <- tuneRF(x = x, y = y,
       stepFactor = 0.9,
       plot = TRUE,
       ntreeTry = 100,
       trace = TRUE,
       improve = 0.05)
## mtry = 5  OOB error = 0.27% 
## Searching left ...
## mtry = 6     OOB error = 0.08% 
## 0.7142857 0.05 
## mtry = 7     OOB error = 0.08% 
## 0 0.05 
## Searching right ...
## mtry = 4     OOB error = 0.58% 
## -6.5 0.05

Aparentemente, 6 é um bom candidato ao valor de \(m_{try}\)

#set.seed(093180)
set.seed(998451)
rf1 <- randomForest(WR_Grp~.,data=train1, 
                    ntree = 100, 
                    mtry = 6, 
                    importance = TRUE,
                    proximity = TRUE)
rf1
## 
## Call:
##  randomForest(formula = WR_Grp ~ ., data = train1, ntree = 100,      mtry = 6, importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 100
## No. of variables tried at each split: 6
## 
##         OOB estimate of  error rate: 18.33%
## Confusion matrix:
##           de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10 class.error
## de 0 a 4         1        7       15         0  0.95652174
## de 4 a 6         0      294      297         0  0.50253807
## de 6 a 8         0      111     1777         4  0.06078224
## de 8 a 10        0        0       43        54  0.44329897
RF_importance1 = randomForest::importance(rf1)[order(randomForest::importance(rf1)[,1], decreasing = TRUE), ]
plot(rf1)
legend('topright', colnames(rf$err.rate), col=1:5, fill=1:5)

RF1 = data.frame(variables = rownames(RF_importance1), importance = RF_importance1[,6])
RF1 = RF1[order(RF1$importance, decreasing = TRUE),]
rownames(RF1) <- NULL

hc6 <- highchart() %>%
  hc_add_series(data = RF1$importance, 
                type = "bar",
                name = "Importância",
                showInLegend = FALSE,
                tooltip = list(valueDecimals = 2, valuePrefix = "", valueSuffix = "")) %>%
  hc_yAxis(title = list(text = "Importância"), 
           allowDecimals = TRUE, max = 200,
           labels = list(format = "{value}")) %>%
  hc_xAxis(title = list(text = "Fatores"),
           categories = RF1$variables,
           tickmarkPlacement = "on",
           opposite = FALSE) %>%
  hc_title(text = "Importância por fator - Random Forest",
           style = list(fontWeight = "bold")) %>% 
  hc_subtitle(text = paste("")) %>%
      hc_tooltip(valueDecimals = 2,
                 pointFormat = "Importância: {point.y}")%>%
                 #pointFormat = "Variável: {point.x} <br> Importância: {point.y}") 
      hc_credits(enabled = TRUE, 
                 text = "Fonte: IMDB/KAGGLE. Elaboração: Ewerson Pimenta.",
                 style = list(fontSize = "10px")) %>%
  hc_exporting(enabled = TRUE, filename = "F4-missing-Pimenta")
#hc <- hc %>% 
#  hc_add_theme(hc_theme_darkunica())
hc6

Predição e matriz de confusão - train data

library(caret)
p1 <- predict(rf1,train1)
head(p1)
##         1         2         3         4         5         6 
##  de 6 a 8  de 6 a 8  de 6 a 8 de 8 a 10  de 6 a 8  de 6 a 8 
## Levels: de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10
head(train1$WR_Grp)
## [1] de 6 a 8  de 6 a 8  de 6 a 8  de 8 a 10 de 6 a 8  de 6 a 8 
## Levels: de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10
confusionMatrix(p1, train1$WR_Grp)
## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10
##   de 0 a 4        23        0        0         0
##   de 4 a 6         0      591        0         0
##   de 6 a 8         0        0     1892         0
##   de 8 a 10        0        0        0        97
## 
## Overall Statistics
##                                      
##                Accuracy : 1          
##                  95% CI : (0.9986, 1)
##     No Information Rate : 0.7269     
##     P-Value [Acc > NIR] : < 2.2e-16  
##                                      
##                   Kappa : 1          
##  Mcnemar's Test P-Value : NA         
## 
## Statistics by Class:
## 
##                      Class: de 0 a 4 Class: de 4 a 6 Class: de 6 a 8
## Sensitivity                 1.000000           1.000          1.0000
## Specificity                 1.000000           1.000          1.0000
## Pos Pred Value              1.000000           1.000          1.0000
## Neg Pred Value              1.000000           1.000          1.0000
## Prevalence                  0.008836           0.227          0.7269
## Detection Rate              0.008836           0.227          0.7269
## Detection Prevalence        0.008836           0.227          0.7269
## Balanced Accuracy           1.000000           1.000          1.0000
##                      Class: de 8 a 10
## Sensitivity                   1.00000
## Specificity                   1.00000
## Pos Pred Value                1.00000
## Neg Pred Value                1.00000
## Prevalence                    0.03726
## Detection Rate                0.03726
## Detection Prevalence          0.03726
## Balanced Accuracy             1.00000

Predição e matriz de confusão - test data

p2 <- predict(rf1,test1)
head(p2)
##         1         2         3         4         5         6 
##  de 6 a 8  de 6 a 8  de 6 a 8 de 8 a 10  de 6 a 8  de 6 a 8 
## Levels: de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10
head(test1$WR_Grp)
## [1] de 6 a 8 de 6 a 8 de 6 a 8 de 6 a 8 de 6 a 8 de 6 a 8
## Levels: de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10
confusionMatrix(p2, test1$WR_Grp)
## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  de 0 a 4 de 4 a 6 de 6 a 8 de 8 a 10
##   de 0 a 4         0        0        0         0
##   de 4 a 6         2      113       41         0
##   de 6 a 8         4      139      753        15
##   de 8 a 10        0        0        5        12
## 
## Overall Statistics
##                                           
##                Accuracy : 0.81            
##                  95% CI : (0.7853, 0.8329)
##     No Information Rate : 0.7371          
##     P-Value [Acc > NIR] : 1.083e-08       
##                                           
##                   Kappa : 0.4519          
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: de 0 a 4 Class: de 4 a 6 Class: de 6 a 8
## Sensitivity                 0.000000          0.4484          0.9424
## Specificity                 1.000000          0.9483          0.4456
## Pos Pred Value                   NaN          0.7244          0.8266
## Neg Pred Value              0.994465          0.8502          0.7341
## Prevalence                  0.005535          0.2325          0.7371
## Detection Rate              0.000000          0.1042          0.6946
## Detection Prevalence        0.000000          0.1439          0.8404
## Balanced Accuracy           0.500000          0.6984          0.6940
##                      Class: de 8 a 10
## Sensitivity                   0.44444
## Specificity                   0.99527
## Pos Pred Value                0.70588
## Neg Pred Value                0.98594
## Prevalence                    0.02491
## Detection Rate                0.01107
## Detection Prevalence          0.01568
## Balanced Accuracy             0.71986

Nº de nós nas árvores

hist(treesize(rf1),
     main = "Nº de nós por ávore",
     col = "green")

Extração de uma única árvore \(Árvore \space n_{tree}=1\)

datatable(getTree(rf1, 1, labelVar = TRUE),  
          options = list(searchin = TRUE, pageLength = 5))
remove(INFO)
## Warning in remove(INFO): objeto 'INFO' não encontrado

\(Árvore \space n_{tree}=100\)

datatable(getTree(rf1, 100, labelVar = TRUE),  
          options = list(searchin = TRUE, pageLength = 5))
remove(INFO)
## Warning in remove(INFO): objeto 'INFO' não encontrado

Gráfico de escala multidimensional da matriz de proximidade

MDSplot(rf1, train1$WR_Grp, pch=20)

Plot do modelo rpart - Recursive Partitioning and Regression Trees - personalizando automaticamente a partir do gráfico para o tipo de resposta do modelo.

rp <- rpart::rpart(formula = WR_Grp~.,data=test1)
#rpart::plotcp(rp)
rpart.plot(rp)

rpart.plot.version1(rp)

Tempo de execução

proc.time() -ptm
##    user  system elapsed 
##   53.91    3.20   60.72